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The dynamics of vortex solitons in a BEG superfluid is studied. A quantum lattice-gas algorithm 
(localization-based quantum computation) is employed to examine the dynamical behavior of vortex 
soliton solutions of the Gross-Pitaevskii equation (^'^ interaction nonlinear Schroedinger equation). 
Quantum turbulence is studied in large grid numerical simulations: Kolmogorov spectrum associated 
with a Richardson energy cascade occurs on large flow scales. At intermediate scales a A:~^ power 
law emerges, in a classical- quantum transition from vortex filament reconnections to Kelvin wave- 
acoustic wave coupling. The spontaneous exchange of intermediate vortex rings is observed. Finally, 
at very small spatial scales a /c"^ power law emerges, characterizing fluid dynamics occurring within 
the scale size of the vortex cores themselves, a characteristic Kelvin wave cascade region. Poincare 
recurrence is studied: in the free non-interacting system, a fast Poincare recurrence occurs for 
regular arrays of line vortices. The recurrence period is used to demarcate dynamics driving the 
nonlinear quantum fluid towards turbulence, since fast recurrence is an approximate symmetry of the 
nonlinear quantum fluid at early times. This class of quantum algorithms is useful for studying BEG 
superfluid dynamics over a broad range of wave numbers, from quantum flow to a pseudo-classical 
inviscid flow regime to a Kolmogorov inert ial subrange. 
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I. INTRODUCTION 

Quantum computation is a newly emerging technique 
within the field of information processing. Feynman's 
conjecture [l], now over a quarter century old, that a 
quantum mechanical computer could efficiently simulate 
physics, principally a many-body quantum system, has 
not yet been physically demonstrated on a quantum com- 
puter. However, quantum algorithms (designed for future 
quantum computers) run today on supercomputers do 
provide a new numerical paradigm for performing stable 
and faithful simulations of quantum systems. Here we 
use a quantum lattice-gas algorithm to model a Bose- 
Einstein condensate (BEG) superfluid in the notoriously 
difficult turbulent regime. The quantum lattice-gas algo- 
rithm is strictly unitary, no dissipative terms are added 
to the reversible dynamics. Thus, the quantum turbu- 
lence observed in the numerical simulation is rigorously 
obtained. 

Someday large-scale quantum simulators should pro- 
vide experimental realization of the Feynman conjecture. 
Until such a time, one can learn much from qubit rep- 
resentations of local dynamics occurring in a quantum 
simulator, particularly from the simplest qubit models 
(perhaps even bootstrapping toward a viable quantum 
computing architecture). In this Letter, we employ a sim- 
ple quantum computational model of a BEG superfluid, a 
minimal model of vortex- vortex {e.g. reconnection) and 
vortex-phonon interactions. Its unitary representation is 
different than the Bose-Hubbard model [2 , which has 
served as the basis of quantum Monte Garlo simulations 
of trapped BEG superfiuids [3 . No low-order Hamilto- 



nian, such as the Bose-Hubbard Hamiltonian accounting 
for particle hopping and particle-particle interactions, is 
specified a priori to generate the dynamics. Instead, ex- 
act dynamics is represented by basic unitary logic gates 
acting at the small lattice cell-size scale. An effective 
Hamiltonian dynamics naturally emerges at large scales. 
Thus the quantum lattice-gas model tested here is a high- 
energy model of the quantum particle dynamics that ef- 
fectively and correctly behaves as a BEG superfluid on 
the large (low-energy) scale. Large grids are employed 
to obtain good numerical predictions, typically 1024^ or 
larger, although accurate numerical results can be ob- 
tained in some cases for grid size an order of magnitude 
smaller. 

Previously, we demonstrated the practicality of our 
quantum lattice-gas computational model by exploring 
instabilities of soliton wave trains in 2+1 dimensions gov- 
erned by a nonlinear Schrodinger wave equation [4] . Gon- 
sidering vortex solitons in BEG superfiuids, we now focus 
on the behavior of vortex lattices as a means to compre- 
hend quantum turbulence in 3+1 dimensions. We use the 
small vortex lattices (quadrupolar arrangements of linear 
vortex solitons) as the basis of initial conditions to study 
the subsequent behavior a BEG superfluid flow. These 
vortex lattices have zero total angular momentum, in 
contradistinction to quantized vortex lattices with fixed 
nonzero angular rotation in one direction that have been 
used as a marker of superfluid flow within a normal fluid 
background [5]. 

The findings reported in this Letter are: 

1. Vortex lattices with multiple orthogonal angular ro- 
tation directions evolve into highly tangled turbu- 
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lent like flows. Energy spectra is consistent with a 
^-5/3 Kolmogorov power law for vortex- vortex flow 
morphology for small wave numbers. 

2. At large wave numbers, the energy spectra is con- 
sistent with a power law. 

3. There is a well-deflned cross-over region with its 
own characteristic power-law behavior. 

4. In free quantum fluids, or BEC superfluids with 
very weak nonlinear coupling, turbulent flow with 
tangled vortex fllaments can quickly emerge but 
then just as quickly untangle. Thus, fast Poicare 
recurrence occurs in a weak BEC superfluid, al- 
beit with remnant Kelvin wave (twisted vortex soli- 
tons). 

5. A highly nonlinear mechanism for vortex- vortex in- 
teraction by exchanging vortex rings occurs in a 
BEC superfluid but also occurs in a free quantum 
fluid (linear Schroedinger wave equation). 

A. Application of measurement-based quantum 
computing 

There are different quantum informational represen- 
tations of quantum computation including the standard 
quantum circuit model and adiabatic quantum comput- 
ing. Type-II quantum lattice gas models use state lo- 
calization to boost the usefulness of quantum algorithms 
for physical simulation [H [3 [8]. Furthermore, type-II 
quantum computational models were the flrst examples 
of "one-way" measurement-based quantum computing 
relying on a parallel array of locally entangled cluster 
states, an antecedent to cluster-based quantum compu- 
tation [QlIIOlIIIlIIl]. 

We employ a representation that relies on 2-qubit state 
localization "self-steering" the outcome of the quantum 
computation. Remarkably this mechanism gives rise to 
quantum fluidic behavior at the large scale, i.e. the non- 
linear interaction term in the Gross Pitaevskii equation, 
the effective equation of motion of a low-temperature 
BEC superfluid. With this representation, previously we 
have predicted solutions to a number of nonlinear classi- 
cal and quantum systems [4[ |13l [T4[ |15l [16] . An advan- 
tage of our approach over the standard computational 
physics CP solvers is that the simple unitary collide- 
stream-rotate operations give rise to an algorithm that 
approaches pseudo-spectral accuracy [TT in much the 
same way as the simple collide-stream steps of a lattice 
Botlzmann algorithm approach pseudo-spectral accuracy 
for fluid dynamics simulations [I8l [19] . 

B. Overview of the modeling approach 

Quantum flelds are discretized on a cubic lattice. Lo- 
cated at each lattice site is a local cluster state, in the 



simplest case just a single pair of qubits. The excited 
state, denoted by logical "1", of a qubit \q{xn,t)) en- 
codes a particle at the lattice site Xn at time t. The local 
ket in the Fock basis of the qubit pair is 

1 

W= E V^g.'k^O=^o|00)+^t|01)+^J10)+^u|ll)- 

q,q'=0 

(1) 

The local Fock states |01) and |10) encode the spin states 
It) and II) and so the arrows as subscripts denote the re- 
spective amplitudes of those states. The mean-fleld ten- 

sor product |^(t)) = (2)n=i l^i^nj t)), where L is the (lin- 
ear) lattice size, represents the state of the quantum sys- 
tem. Furthermore, to analytically recover the CP equa- 
tion it is suflicient to consider just the one-body sector, 
hence the 2-component fleld 

is computed on the lattice. The reduction from ([T]) to 
([2| is not an approximation. Q is an exact representa- 
tion of the quantum dynamics so long as the simulation 
is carried out in the one-body sector. Finally, the BEC 
wave function (j) is determined as the sum of the ampli- 
tudes (j) = ip^ This, of course, retains the important 
effect of quantum interference. 

The evolution of is determined by j-j qubit-qubit 
interactions (collide), free motion of the amplitudes along 
the cubic lattice (stream), and qubit phase shifts to 
model the well known "Mexican hat" interaction po- 
tential (rotation). The qubit-qubit collisions are gener- 
ated by a collision Hamiltonian = Jctx^ where J is 
the |-| coupling energy, ax a Pauli spin operator, and 
the collision time r corresponds to the quantum logic 
gate time, chosen such that = f • Free streaming of 
qubit states on the lattice (emulating the motion of par- 
ticles in space) is generated by the stream Hamiltonian 
Hs = — ^^Zliattice ■ ^ acting on the qubits. c is a 
streaming velocity along lattice directions. 

Exploiting the fact that [i^s,^c] 7^ 0^ we use an inter- 
leaved compositional product of 2-qubit quantum gates 
generated and He whereby an effective Hamiltonian 
emerges in the scaling limit modeling the non-relativistic 
free particle Hamiltonian Hq. Formally denoting this 
compositional product by the symbol o, the Hamiltonian 
Ho = Tt[Hs o He] + • • • 7 where in our model 

m = ^ in lattice units and where the right arrow denotes 
a scaling-limit mapping. Thus Ho generates the evolu- 
tion of a free scalar quantum fleld. If this does not seem 
clear at the moment, do not fret because an explanation 
is presented in the section following immediately below. 

Phase rotation, inducing nonlinear particle-particle in- 
teractions, is generated by i^int(I^P) = (^|(/>P — l) |(/>P, 
where the quantum logic gate r' is chosen such that 
^ ^ 1. In the low energy limit, the local evolution 
is effectively 

^{t + 1, x) = e-^^°"/^ e-^^^-^"'/^ ^{t, x) + • • • , (3) 
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modeling (25) in the lowest order fluctuations. A deriva- 
tion of Hqp ^ Ho -\- i^int presented in [20 follows from 
the typical scaling argument used in kinetic lattice gases, 
which we also review below. 



II. QUANTUM LATTICE-GAS ALGORITHM 



Conservative two-qubit quantum gates have the form 



m 



conservativ 



a 
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(4) 



where the components in the zero-quantum subspace are 



a member of SU(2), 



VIZ. 



C V 



e SU(2). Only the 



local 2-component fleld's complex amplitudes ip^{x) and 
i/jl (x) are quantum mechanically entangled by the action 
of Q. Particle motion and particle-particle interactions 
are faithfully emulated strictly using quantum logic gates 
of the form of Q. Furthermore, to describe the quantum 
lattice-gas algorithm, it is sufficient to consider only the 
single-particle sector of the full quantum Hilbert space. 
In this way, the algorithmic treatment becomes straight- 
forward to describe using only 2x2 matrices that repre- 
sent the SU(2) subspace of Q. 

The justiflcation for this reduction is given in Ref. [13]. 
The quantum gate dynamics conserves particle number 
and consequently the effective H^p in Q commutes with 
the particle number operator. Thus, H^p is block diag- 
onal over the n-body sectors of the Hilbert space. One 
is justified to run a simulation in any one of the n-body 
sector, for < n < 2L^, which is an exact representation 
of all the relevant quantum dynamics in that sector, in- 
dependent of the dynamics in all the other sectors. The 
advantage of limiting the simulation to the one-body sec- 
tor is that the algorithmic complexity scales linearly with 
the number of qubits and the one-body sector is sufficient 
to capture all the relevant physics in the BEC. It is for 
this reason that the type-II quantum algorithm can be 
implemented on a classical supercomputer while retain- 
ing its exactness as a quantum mechanical simulation. 

We now describe the quantum algorithm by dealing 

with the zero-quantum subspace part of the 

quantum logic gate operators Q. The SU(2) quantum 
operator C = e*^^^*^^"^^) that acts locally at every point 
X by the map 

local collision : ^/^(x) Cil){x). (5) 

The well known Pauli matrices are cjx — {^i o) ' 

(jy = ^ = -ij ' The complex scalar density 

(j) = {1^1) • = tp^ -\- ipi is conserved by ([5|, and con- 
sequently the probability is also conserved locally. 



Local equilibrium {i.e. ip = C ip) occurs when the am- 
plitudes are equal {ip^ = i^i)^ but in general such a local 
equilibrium is then broken if a spinor component is dis- 
placed in space by the vectorial amount Ax. To conserve 
probability, we admit only complementary displacements 
of the field components, induced by the stream operators 
of the form 



- e^^^^ n. 



(6a) 



where n = ^(1 — cFz) and n = i(l + ctz)- Although 



the application of (6a) usually breaks local equilibrium 
induced by ([5|, with the appropriate boundary condi- 
tions, for example periodic boundary conditions, ( l6a| ) is 
guaranteed to conserve the total density J d^x ip{xj~8ind 
in turn the total probability J (fix [\ip^{x)\'^ + \ipi{x)\'^) . 
To construct a quantum algorithm using a combination 
of the operators C and SAx,-f^ for 7 = or 1, and the re- 
spective adjoints, we restrict our considerations to those 
combinations which are close to the identity operator. 
Our basic approach uses the interleaved operator 



(7) 



as the basic building block of the quantum algorithm. For 
example, an evolution operator for the 7th component is 



(8) 



where ^ ^ where N is the grid resolution {i.e. N is the 
number of grid points along one edge of the simulation 
volume). (|8| represents the three aspects of a type-II 
quantum algorithm: stream, collide, and state reduction. 
In dimensionless units (c = 1), note that ~ Ax^ ~ 
At. This evolution operator is spatially dependent only 
through local state reduction Q: 

i){x,t^At) = u^[n]iij{x,t). (9) 

(|9| specifies the nonlinear quantum lattice gas model. 

A. Effective field theory 

The R.H.S. of (|9| can be expressed as a finite-difference 
of t/j. Taylor expanding the R.H.S. in e, one obtains the 
following quantum lattice gas equation 



+ At) tp{x,t) 



le 



1 



+ ^ 



ipix, t) 



K + (T,)V3^(f,t) + 0(e4), 



(10) 



where 7 = or 1. Since the order £ error term in (10) 



changes sign with 7, we can induce a cancelation of this 
error term using a symmetric evolution operator 



Uo 



(11) 
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that is invariant under field component interchange. 
Therefore, the fohowing quantum map t + At) = 
U[Q{x)]ilj{x^t) leads to the quantum lattice gas equation 



1 



^{x,t)^0{e^). 

(12) 

Now, in the low energy scaling limit, we have 
^ , At)- ipix, t)] dtipix, t). Therefore, divid- 
ing both sides of (12) by e^, the quantum lattice gas 
equation is 



(13) 



in the low energy and low momentum limit. Finally, since 
(/) = a-\- taking the density moment gives the effective 
scalar field equation 



idt^ = -\/^^^n^^O{e^), 



(14) 



which is the Schroedinger wave equation with m = ^ for 
= 1, so long as |Axp = At = £. From the order of the 
error term in ([l4| , the Taylor expansion predicts that the 
quantum algorithm is second order convergent in space. 

The low energy effective Hamiltonian that is the model 
generator of the evolution, U = e^^^^^ff/^^ follow- 
ing 



^eff 



-^\/^^hQ{x) 
2m 



0{At,Ax^), (15) 



where we have written the quantum diffusion coefficient 
as 4^ = — . This is the nonlinear GP Hamiltonian since 

At m 

hVt{x) = g\(j){x)\^ — 1, where g the on-site interaction 
energy. 



B. Matching BEC experiments 

BEG optical lattice setups tuned to a Feshbach res- 
onance can simulate either repulsive or attractive g = 
with tto is the scattering length [2l]. In the simu- 
lations presented below a repulsive nonlinear interaction 
is used. We have been careful not to choose g too large, 
as (25) is appropriate for weak scattering; a large g trig- 



gers the occurrence of three-body molecular interactions, 
violating (25) and destroying the BEG. In setups with 
^ 10^ atoms, g < 100 is accessible. So a quantum lattice 
gas simulation is realistic so long as g < 100. 



We use a spatial rescaling parameter ^/a so that the vor- 
tex core can be resolved on the computational lattice 
with sufficient numerical accuracy. (16) is modeled by 
using = — 1 in (11), and this corresponds to an 
interaction Lagrangian density £int = a^*^(l — g(j)*(j)). 
Then (ITsl) becomes the nonlinear Hamiltonian for the 



GP equation (16). With the appropriate choice of non- 
linear coupling g and normalization of the wave function, 
physically the energy of a constant external potential can 
cancel the internal interaction energy, at least in the re- 
gion of bulk flow. This is the type of solution that I will 
explain here. 

A solution for the wave function of the vortex soliton 
is found by separation of variables in polar coordinates. 
Inserting (j){r^6) = /(r)e*^^, for integer n, into (16) then 
gives 



d?f{r) , 1 df{r) r? 



adr^ 



ar dr 



.f{r)+{l-gf{rf) f{r) = 0, 

(17) 



which can be solved for any winding number. For the 
simplest n = 1 and g = I case, a Pade approximant of 
the spatially scalable form is 



fir) 



Then, rescaling ^/a^ 



' llar2(12 + 


ar'^) 


384 + ar2(128 


+ llar2) 


^ r, we have 





(18) 



R"{r) + -R'{r) - \R{r) + (l - R{rf) R{r) = 0, (19) 

where R{r) = / ^-^^ is the solution of the radial part 

of ( [25| ) found by Berloff [22]. Rescaling allows one to 
arbitrarily resolve the vortex core on the computational 
lattice to achieve sufficient numerical accuracy. 

([l6| satisfies the free field condition in the bulk 
when ^ (far away from any vortex singularity). 

This is called hulk point normalization. That is, since 
lim^^oo /(^) = 1, to satisfy the free field condition we 
must choose g = 1. Thus with y^r ^ r, at the end of 
the day the equation we are going to numerically solve is 



zS,(/) = -V'(/)+(|(/)p-l)(/) = 



(20) 



if we were to model a single vortex soliton. Unfortu- 
nately, in a box with periodic boundary conditions we 
cannot make do with just a single vortex soliton. We 
need at least four of them. 



III. VORTEX SOLITONS 

A. Single line soliton 

We first seek a line vortex steady state solution of the 
rescaled GP equation 



1 



V^(/>+(^|(/>p-l)(/> = 0. 



(16) 



B. Arbitrary coupling strength via wave function 
normalization 



The GP equation is invariant under arbitrary wave 
function normalization provided one also rescales the 
nonlinear interaction: 



idt^ = -V^^^a{g\^\^ -l)ct) 



(21) 
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• and g 



is invariant under (j) {aa^ Li)~ 
where a is the spatial scahng parameter in (18) and N is 
the number of vortex lines, e.g. A/" = 4 in (22), L is the 

grid size, and the factor a = a~^L~i (J d^x\(j){x)\'^) ^ ^ 
1 accounts for the excluded volume due to the vortex 



C. Quadrupolar configuration (4 line solitons) 




FIG. 1: A slice at z = Zo of the magnitude p(x, y, Zo) (top, 
upside down) and phase 0(x,y,Zo) (middle) and phase contours 
(bottom) of the wave function for a vortex quadrupole, the product 
of 4 vortex soliton solutions on a grid of size L = 160. The density 
p(x,y) = |0(x, 2/)p and so -sfp^x^/y^Zo) 1 away from a vortex core 
(in the bulk). From the phase diagram, plotted — tt < 6(x,y,Zo) < 
TT, going around any contour in the z = Zo plane that encloses 
a single vortex singularity accumulates a phase of ±27r radians. 
With N = 4 line vortices one can accommodate periodic boundary 
conditions in the phase. 



A simple initial condition that ensures periodicity is 
four symmetrically displaced vortex line solitons (parallel 



to the z-axis for the time being) in product form 

cj,{x,y) = /(r++)e*^++ X /(r+_)e-*^+- 

X f{r_+)e-''-+ X /(r__)e*^-- (22a) 

= /(r++)/(r+_)/(r_+)/(r__) 

X e'^'^++-'+—'^-++'--\ (22b) 

where the radial distance from a vortex line along the 
2;-axis is 



r++(x,y) = {x-Xo 

r^-{x,y) = {x - Xo 

r-^{x,y) = {x-Xo 

r — {x,y) = {x-Xo 



S)^^{y-yo^S)^ (23a) 

S)'^{y-yo-Sf (23b) 

Sf^iy-yo^Sf (23c) 

S)'^{y-yo-S)^ (23d) 



The size of the vortex quadrupole is |2^|, sign(5) = ±1 
is called the polarity of the vortex quadrupole, and its 
center is {xo^yo)- The phase angles are 
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(24a) 
(24b) 
(24c) 
(24d) 



The magnitude and phase of (22) are plotted in Fig. [T] 
with S = ^ and a = 0.1 and L = 160, demonstrating the 



periodicity of (22). We shall use such quadrupole vortex 
lattice configurations aligned along orthogonal principal 
lattice direction to represent initial conditions for numer- 
ical simulations. 



IV. QUANTUM TURBULENCE 

Turbulence remains one of the greatest unsolved prob- 
lems of classical physics, even in the incompressible limit. 
The complexities of classical turbulence are further com- 
pounded by the inability to accurately define a classical 
"eddy." This creates ambiguity in the phenomenology of 
the inertial cascade [23 and the self-similarity in the dy- 
namics of the smaller eddies that leads to the famous 
Kolmogorov k~^^^ energy spectrum p4| [25 ] . 

Feynman proposed that the superfiuid turbulent state 
consists of a tangle of quantized vortices [26]. Quantum 
turbulent studies, both experimental [27 and computa- 
tional [28j ^29^ , are being pursued vigorously-not only for 
their intrinsic importance but also in the hope of shed- 
ding light on classical turbulence. This is consistent with 
behavior observed in our quantum algorithms simulations 
of a BEC. In contrast to classical eddies, a quantized vor- 
tex is a well defined stable topological line defect in the 
phase of the complex scalar quantum wave function and 
a soliton in the magnitude of the wave function. The vor- 
ticity is generated in a multiply-connected region where 
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FIG. 2: Starting with = 12 vortex lines on a 512^ lattice. 
Vortex tubes at t = 200K (top) show an onset of a Kelvin wave 
instability. Tangled vortices are observed, even when i^int ~ 0, at 
t = 3.3K (bottom). Remarkably, one observes many vortex rings 
mediating the vortex line-line interactions. 

the circulation is quantized: the phase of the complex 
scalar wave function along a contour is accumulated in 
integer multiples of 27r for any contour enclosing a quan- 
tum vortex center. Early superfluid turbulence studies 
focused on a regime describable by Landau's two-fluid 
equations that modeled the interaction between the in- 
viscid superfluid and a viscous normal fluid [30^ and un- 
derstood in terms of the motion of quantized vortices due 
to the Magnus force In constrast, a BEG superfluid 
comprises a macroscopic number of integer spin particles 
in effectively a zero temperture ground state, and its dy- 
namics below the critical temperature is remarkably well 



described by the Gross-Pitaevskii (GP) equation [32j[33], 
which in normalized form is 

idt<P=-V^<l>+M^ -1)<P, (25) 

where (j) is the mean-field BEG wave function. It is well 
known |34| that the Madelung transformation (j) = y/p e*^ 
on the GP equation results in compressible inviscid fluid 
equations for the density p = and velocity v = 2V6>, 
with the appearance of quantum pressure terms in the 
momentum and energy equations. 



A. Viewing quantum turbulence 

As mentioned above, products of quadrupolar line soli- 
tons depicted in Fig. [l] may be shifted and rotated about 
any direction to build up unstable yet periodic initial con- 
ditions. Here we report on a simulation initialized with 
3 orthogonal sets of 4 vortex centers (total of = 12 
vortex lines) on a 512^ and 1024^ grids. For A = 12 
superfluid simulation shown in Fig. |2j transverse vortex 
waves (helical perturbations of the vortex tube away from 
its cylindrical shape) known as Kelvin waves onset on the 
filamentary core hy t = 200 time steps and subsequent 
become unstable, the initial drive force of the energy cas- 
cade [35 (large waves coupling to smaller waves). There 
is a rapid tangling of the quantized vortices (by t ~ 3K) . 

One might expect a stark difference between the ob- 
served quantum turbulence arising from the compressible 
inviscid fluid equations of the GP equation without vis- 
cous shear dissipation, since unitary quantum dynamics 
is Hamiltonian, and classical turbulence arising from the 
incompressible Navier- Stokes equation with viscous shear 
dissipation, since viscous hydrodynamics system is non- 
Hamiltonian. There are clear similarities for small wave 
numbers. 



B. Spectrum for incompressible kinetic energy 

The approach we follow is to decompose the conserved 
energy into its kinetic, quantum and internal energy parts 

as follows Etot = ^kin + Ec^u + ^int = COUSt. 

(26) 

Vortex reconnections through Kelvin wave instabilities 
contribute to the Richardson cascade [23 (large vor- 
tices break into smaller ones) leading to the famous Kol- 
mogorov spectrum. The presently held notion is that vor- 
tex filaments reduce down into smaller and smaller loops 
and the high-/c oscillations of the vortex center couple 
to elementary phonon excitations at the healing length 
giving rise to a quantum dissipation scale. Vortices are 
destroyed by the emission of sound waves (compressible 
excitations) that escape across the space. 
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spatial cut-offs, 8.44 < r < 19.71. The A:"^-^^ power law 
occurs on spatial scales larger than vortex core size. For 
very large wave numbers k > 300 (with a cut-off of about 
just 2 lattice cell sizes), the spectrum all together drops 
off the chart, as expected. 



V. POINCARE RECURRENCE 




FIG. 4: In N = 12 simulation, untangled vortices are observed 
t = 21K (bottom). The inital state recurs after a turbulent state. 
An ordered state at t = 21K deterministically returns to the initial 
state untangling turbulence, a cascade that cycles at intervals of 
tp = 21K on a 512^ lattice. 



FIG. 3: Power spectrum of the quantum fluid's incompressible 
kinetic energy (top). There are three regions characterized by dif- 
fering power laws displayed on the vortex soliton spatial proflle 
(bottom). Numerical data (dots) is from a supercomputer simula- 
tion of a quantum lattice gas on a 1024^ grid. Kolmogorov (black), 
transition (green), and core interior (red) regimes are shown. 



A Kolmogorov spectrum is observed in incompressible 
and quantum energies of the BEG superfluid for /c < 30, 
see Fig. [3j The emergence of full developed quantum 
turbulence is plotted at time t = 27. 7 K for a 1024^ lat- 
tice starting with = 12 vortex solitons (1 quadrupole 
per spatial direction). The measured power law k~^-^^ 
(black) may suggest that the theoretical Komolgorov 
power law describes the spectrum for /c < 20, although 
the fit is not excellent. A new power law k~^-^^ (green) 
emerges in a transition region from 30 < /c < 70. After 
/c > 70, the observed power law k~^-^^ (red) agrees with 
the theoretical prediction of k~^ and excellently fits the 
data in this region. The power law fits were all computed 
using linear regression. 

The bottom plot in Fig. [3] shows the three power 
regimes in a spatial view by overlaying the cut-off lengths 
on the vortex core profile. The k~^ power law character- 
izes fluid dynamics within the vortex core itself, with an 
upper cut-off scale measured to occur at r = 8.44 (in 
units of lattice sites). The k~^-^^ power law occurs near 
the boundary of the vortex core, with lower and upper 



In Hamiltonian systems, the dynamics must be invert- 
ible, so it is possible to observe a reverse cascade, and on 
time scales shorter than otherwise expected; vortex tubes 
untangle and reform by the absorption of sound waves to 
recover a configuration close to a configuration that oc- 
curred earlier in time in the flow. Recursion arising from 
the Hamiltonian system is observed in animations of the 
flow. This occurs in the limit of vanishing nonlinear in- 
teraction, ^ ~ 0, the vortex solitons completely untangle, 
as evidenced in Fig. |4] (by t = 21i^), when the internal 
energy in (26) satisfles £^int <^ £^kin, ^qu- Surprisingly in 



3+1 dimensions, the Poincare recursion time for the GP 



equation (25) can be extremely short. Furthermore, fast 



Poincare recurrence occurs only for simple initial condi- 
tions. 

In a strictly unitary quantum algorithm, it is possi- 
ble to observe recurrence because the numerical method 
is information preserving. In a simulation in a flnite- 
sized box with periodic boundary conditions, when kelvin 
waves excite phonons (incompressible energy changing to 
compressible energy) this is indeed acts like a dissipative 
mechanism on the large scale and it gives rise to the 
characteristic k~^ power law that we observe at large k. 
However, the excited phonons cannot escape off to infln- 
ity since the simulation is in a flnite-sized box. At the 
point half way through one Poincare recurrence cycle, all 
the available phonon states are essentially fllled and the 
quantum fluid can no longer transfer incompressible en- 
ergy into these compressible modes. It is at this point 
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that the net transfer of energy reverses, and an inverse 
cascade is observed. Although, we refer to the Poincare 
recurrence time as fast, this is a relative term. On a small 
512^ grid, one would have to run through over ~ 20, 000 
updates to see recurrence and this is a rather long simula- 
tion time. Previously the computational resources where 
not available for do this kind of long run time, particular 
for pseudo-spectral algorithms that slowdown on large 
grids. Furthermore, using complicated initial condition, 
such as the Taylor-Green profile [34 , does not admit the 
shortest Poincare cycle time in the first place. Finally, 
as far as we can tell, tests of turbulent superflow in the 
case with g very small have not been conducted. This is 
why recurrence has not been previously observed. 




0,0 

20000 40000 60000 80000 100000 

t (time step) 
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10-^ 



10000 20000 30000 40000 50000 

V (frequency bins) 

FIG. 5: Top: The time evolution of the kinetic energy (blue), 
quantum energy (red), internal energy (green), and total energy 
(black line). The total energy is conserved. A recurrence of tp = 
20. 9X time steps is determined from the envelop of the forward 
energy cascade (dashed). Bottom: Power spectrum indicates all 
frequencies are represented as \E{h>)\^ ~ [z^(z^max — ^)]~^ (red). 



To quantify this phenomenon, one can plot a time his- 
tory of the kinetic, quantum and internal energy parts. 
In the = 12 simulation with i^intd^P) ^ reverse 
cascade (absorption of sound waves) clearly repeats, cy- 
cling at tp ~ 21i^, the Poincare recursion time for grid 
size L = 512. This recurrence time is clearly evidenced 
in the time evolution of the (rescaled) kinetic £^kin and 
quantum Ec^^ energies plotted in Fig. |5] 



The Poincare recurrence theorem states that for 
Hamiltonian systems the solution trajectory passes ar- 
bitrarily close to the initial state provided the evolution 
is followed for a sufficiently long time. While for certain 
maps in two spatial dimensions, like the Arnold Cat Map, 
the Poincare recurrence time can be short, for nearly all 
Hamiltonian systems the recurrence time is so long as 
to be effectively infinite. There have been some analyti- 
cal hints that the NLS equation in 1+1 dimensions could 
have a fast Poincaree recurrence time [36 -but this re- 
sult was not expected to hold in three spatial dimensions. 
From a series of quantum simulations we see the Poincare 
recursion time scales as , where L is the (linear) grid 
size. This arises from the inherent diffusive ordering of 
fluctuations, —St^ 5x^ . 

VI. ISOLATING KELVIN WAVES 

We examined hydrodynamic-scale breaking of the 
quantum-scale time-reversal symmetry of the free sys- 
tem, i^int breaks the Poincare recursion and has a promi- 
nent effect on the dynamics of the quantized line vortices. 
Yet, it remains useful to chart the pathway to turbulent 
configurations at intervals demarcated by the Poincare 
recurrence time of the free system. In the interacting 
GP limit, the fast recursion is broken by nonlinear twist- 
ing (Kelvin waves) riding on the originally linear vortices, 
i^int successively twisting the filamentary centers every 
recursion period; see Fig. [6j The linear vortex tubes 
become tangled but at the free recurrence period they 
do not return to their original linear configuration. In- 
stead, they become twisted and this twisting increases 
with more and more free Poincare cycles. 

At very large times, the BEC manifests pure quantum 
turbulence, characteristic of nonlinear fluid behavior. In 
this numerical simulation, L = 160 and the smallness in 
the nonlinear interaction in ^ is set to ^ = 0.1. To 
sufficiently resolve the vortex core, the scale factor in the 
Pade approximant is set to a = 0.05. A convention of 
unity normalization is used (J \(j){x)\'^dx^ = 1). Poincare 
recurrence in the g ^ limit occurs at tp c:^ 2020i^, 
and this time period is used to sample the wave function 
configurations of the GP quantum system with g ^ 5. 

VII. CONCLUSION 

A quantum lattice-gas algorithm for modeling 
hydrodynamic-scale BEC superfluid flow was presented. 
The method accurately captures the dynamical behavior 
the superfluid even in the difficult case of fully developed 
quantum turbulence. Interacting vortex soliton lead to 
fluid instabilities in the unitary quantum fluid that causes 
an energy cascade in the regime of small wave numbers, 
leading to power law behavior k~3 Kolmogorov turbu- 
lence in classical Navier- Stokes fluids. However, at larger 
wave numbers (approaching the scale of the vortex cores). 







FIG. 6: Kelvin waves seen as twisting of (N = 8) vortex filaments when //int(|</>P) = ~ W^)- the g ^ limit (non-interacting 

particles), there is a fast Poincare recurrence time of tp = 2020. For 5^ ~ 5, vortices at the first few Poincare cycles (t = 0, tp, 2tp, 3tp) are 
plotted (top to bottom). The Kelvin wave twisting in the vortices eventually completely breaks the fast Poincare recurrence. The highly 
tangled vortices, similar to that of Fig. 2, occurring between the Poincare periods are not shown. The asymmetry in the time is due to 
the broken symmetry of the initial condition. 



the spectrum of incompressible kinetic energy transitions 
to an alternate power law unique to quantum fluid flow, 
than finally to a for very large wave numbers (corre- 
sponding to fluid dynamics internal to the vortex cores). 
This power-law behavior has recently been observed in 
experimental observations of superfluid Helium II [37] . 

For weakly nonlinear BEG superfluids, it is possible to 
observe fast Poincare recurrence where at each recurrence 
time, Kelvin waves are observed to emerge with greater 
amplitude at each recurrence time, viz., vortex solitons 
that are successively twisted. 

The quantum algorithm tested here hopefully will be 
implemented on a large quantum computer someday, if 
that device becomes available. At that time, the many- 
body quantum simulation problem can be directly ad- 
dressed numerically using the quantum complexity in- 
herent to the quantum computer. This would include in- 
terference effects across the n-body sectors of the Hilbert 
space caused by the nonlinear particle-particle interac- 



tions, whereas in our present supercomputer-based sim- 
ulation we have restricted ourselves to modeling a one- 
component BEG superfluid in the one-body sector. 
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